%Structual estimation (main file)


clear
tic

choice=1:13;
beta=0.95;
T=100;
length=0:T;
df=beta.^length;
%AI_fs=[51:63,65:66,68:72,76:94,96:119,121:145,147:150];


start=1;
page=start:(start+29);
narray=numel(page);
%narray=s;
%coef=zeros(narray,2);
%ste=zeros(narray,2);
%fv=zeros(narray,1);

fs0=dlmread('fs1.csv');
fs0=sortrows(fs0,[1,2,5]);
ddom0=dlmread('mks1.csv');
ddom0=sortrows(ddom0,[1,2,3]);
FSTA=zeros(size(fs0,1),size(fs0,2),narray);
DOM=zeros(size(ddom0,1),size(ddom0,2),narray);

% Construct other inputs from "regvar*"
bd=fs0(:,3,1);
%nobd=bd(:,ones(T+1,1));
%bdtype=(bd<=25&bd>0)+(bd<=80&bd>25)*2+(bd<=193&bd>80)*3+(bd>193)*4; %bdtot: 25%: 25; 50%: 80; 75%: 193
tmpbd=unique(bd);
tmpbd(tmpbd==0)=[];
bed25=prctile(tmpbd,25);
bed50=median(tmpbd);
bed75=prctile(tmpbd,75);
bdtype=(bd<=bed25&bd>0)+(bd<=bed50&bd>bed25)*2+(bd<=bed75&bd>bed50)*3+(bd>bed75)*4;
%bdtype=[(bd<=bed25&bd>0), (bd<=bed50&bd>bed25), (bd<=bed75&bd>bed50), (bd>bed75)];



textFileName26=sprintf('regvar6_sing.raw');  %the input file  
textFileName27=sprintf('regvar7_sing.raw');  %the input file     
textFileName28=sprintf('regvar8_sing.raw');  %the input file     
textFileName29=sprintf('regvar9_sing.raw');  %the input file    

char6=dlmread(textFileName26); %the matrix of characteristics of stand-alone hospitals: ahaid year bdtot prechs adj_vnd hsa
char7=dlmread(textFileName27);    
char8=dlmread(textFileName28);
char9=dlmread(textFileName29);

allchar=[char6;char7;char8;char9];
allchar=sortrows(allchar,[1,2]); 
pre=allchar(:,4);
act=allchar(:,5);  %extract the current action to compute the probability of the actual choice

% Define cuboid for fminsearch
df_cub=repmat(df,size(FSTA,1),1,narray);
bd_cub=repmat(bd,1,T+1,narray);
bdtype_cub=repmat(bdtype,1,T+1,narray);
gg_cub=zeros(size(DOM,1),T+1,narray);
fchs_cub=zeros(size(FSTA,1),T+1,narray);
ind_cub=zeros(size(FSTA,1),T+1,narray);
idx_sw_cub=zeros(size(FSTA,1),T+1,narray);
eerr_cub=zeros(size(FSTA,1),T+1,narray);


%### UPDATED on 07/04/2020: create variables for market (un)observables
%textFileName_mktvar=sprintf('mktvar.raw');  %the input file 
textFileName_mktvar=sprintf('mktvar_pophhi.raw');  %the input file
mktvar=dlmread(textFileName_mktvar); %the matrix: ahaid year hhi mktcat totpop65
mktvar=sortrows(mktvar,[1,2]);
mktcat_tmp=zeros(size(fs0,1),1);
hhi_tmp=zeros(size(fs0,1),1);
pop65_tmp=zeros(size(fs0,1),1);
for i=1:size(mktvar,1)
    hhi_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,3);
    mktcat_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,4);
    pop65_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,5);
end
hhi_tmp(fs0(:,3)==0)=0;
mktcat_tmp(fs0(:,3)==0)=0;
pop65_tmp(fs0(:,3)==0)=0;
%hhi=hhi_tmp(:,ones(101,1));
%mktcat=mktcat_tmp(:,ones(101,1));
hhi_cub=repmat(hhi_tmp,1,T+1,narray);
mktcat_cub=repmat(mktcat_tmp,1,T+1,narray);
pop65_cub=repmat(pop65_tmp,1,T+1,narray);


eer0=dlmread('eerr1.csv'); %matrix: ahaid year option_given pr1 pr2 ... pr100
EERR=zeros(size(eer0,1),size(eer0,2)+1,narray);   %#### UPDATED on 05/05/2017: (Take log!!!) the output for the expected logit error: ahaid, option given, prob@t=1, prob@t=2, ..., prob@t=100
clear hhi_tmp mktcat_tmp pop65_tmp
clear fs0 ddom0 eer0
for c=1:narray
    s=page(c);
    textFileName1=sprintf('fs%d.csv',s);
    textFileName2=sprintf('mks%d.csv',s);   
    textFileName3=sprintf('eerr%d.csv',s);
    tempEERR=dlmread(textFileName3);
    tempEERR(tempEERR==0)=1;
    EERR(:,:,c)=[tempEERR(:,1:3),zeros(size(tempEERR,1),1),log(tempEERR(:,4:end))]; %#### UPDATED on 03/13/2017:include one zero column before all the non-zero expected errors because the first one flow has no expected error.
    EERR(:,:,c)=sortrows(EERR(:,:,c),[1,2,3]);    
    FSTA(:,:,c)=dlmread(textFileName1);  %the matrix: ahaid year bdtot prechs action_given fs1 fs2...fs100
    DOM(:,:,c)=dlmread(textFileName2);  %the matrix: ahaid year action_given ddom_pre ddom0 ddom1 ddom2...ddom99   
    FSTA(:,:,c)=sortrows(FSTA(:,:,c),[1,2,5]);
    DOM(:,:,c)=sortrows(DOM(:,:,c),[1,2,3]);
    
    %### UPDATED on 06/22/2020: try contructing cubic so as to remove the for loops in fminsearch 
    % ev of mktdom
    gg_cub(:,:,c)=DOM(:,4:end,c);
    % cost-saving per bed
    fchs_cub(:,:,c)=FSTA(:,4:end-1,c)+11;    
    % sunk costs
    fir=FSTA(:,4:end-1,c);
    lat=FSTA(:,5:end,c);
    ind_cub(:,:,c)=(fir~=lat).*lat;
    % switching cost
    idx_sw_cub(:,:,c)=(fir~=lat).*(fir>1);
    % the expected errors
    eerr_cub(:,:,c)=EERR(:,4:end,c);    
end

%hbed=FSTA(:,3,1); %bdtot: 25%: 25; 50%: 80; 75%: 193
%bdtype=(hbed<=25&hbed>0)+(hbed<=80&hbed>25)*2+(hbed<=193&hbed>80)*3+(hbed>193)*4;



%%%%%%%%%%%%%%%%%%%%%%%%%%
%with 26 parameters
%theta0=[-6.752912 -6.311207 -3.715860 -4.035594 -6.650967 -6.566256 -7.581982 -5.019611 -5.121039 -5.779016 -4.164671 -4.502513 0.000369 0.000534 0.000058 -0.000058 0.000528 0.000580 0.000505 0.000105 0.000477 0.000506 0.000357 0.000294 0.062659 -1.8];
%theta0=[-5.869513 -5.330296 -2.948035 -3.377029 -5.751234 -5.651111 -6.713429 -4.106970 -4.216721 -4.870577 -3.188701 -3.664923 0.000197 0.000399 -0.000042 -0.000093 0.000394 0.000460 0.000374 -0.000148 0.000345 0.000378 0.000253 0.000127 0.028334 -1.428617];
 
%{
 theta=[-5.6727427	-5.4148654	-2.8961247	-3.2231362	-5.8298762	-5.7155734 ...
  -6.9408619	-4.1749532	-4.2832556	-4.8858257	-3.2607835	-3.4358166 ...
  0.00195965	0.00363447	-0.00157124	-0.00305057	0.00370564	0.00420194 ...
  0.00413653	-0.001424	0.00309183	0.00308829	0.00197276	0.00138111 ...
  0 0 0 ...
  0.065  0.055  0.045  0.035  -1.6  -1.5  -1.4  -1.3];
%}

%{
 theta=[-5.6727427	-5.4148654	-2.8961247	-3.2231362	-5.8298762	-5.7155734 ...
 -6.9408619	-4.1749532	-4.2832556	-4.8858257	-3.2607835	-3.4358166 ...
  0.00195965	0.00363447	-0.00157124	-0.00305057	0.00370564	0.00420194 ...
  0.00413653	-0.001424	0.00309183	0.00308829	0.00197276	0.00138111 ...
  0 0 0 0 0 0 ...  %### UPDATED on 07/04/2020: HHI for each vendor
  0 0 0 0 0 0 ...
   0 0 0 0 0 0 ...  %### UPDATED on 07/04/2020: totpop65 for each vendor
  0 0 0 0 0 0 ...
   0 0 0 ...  %### UPDATED on 04/23/2020: mktcat1(<=25%) mktcat2(<=50%) mktcat3(<=75%) HHI
  0.065  0.055  0.045  0.035  -1.6  -1.5  -1.4  -1.3];  %### UPDATED on 08/11/2018: the last eight: mktdom (4-1) and swc (1-4)
%}

 
theta0=[-5.6599851	-5.5365167	-3.0538648	-3.0282875	-5.3222437	-5.8405458 ...	
-4.9060924	-3.975923	-4.576892	-4.4572401	-3.6615842	-2.7422904 ...		
0.001533366	0.001852315	0.001090066	0.000577794	0.001817492	0.001923451 ...		
0.001859473	0.000753108	0.001821671	0.001731935	0.001634341	0.001054977 ...		
0.02140	0.17140	0.13873	0.06973	-0.10627	0.00651 ...	
-0.04973	0.02553	0.17727	0.11473	0.15127	0.16700 ...	
0.00639	0.03547	0.01493	-0.00753	0.00508	-0.00014 ...	
-0.02300	0.00213	0.01807	0.02707	0.02067	0.01953 ...	
0.012 0.024 0.6 ...
0.1507247	0.17716169	0.11699648	0.01554261	-2.2107251	-2.4743485	-2.0091092	-1.478018]; %### UPDATED on 05/11/2020: try new starting point



%{
%theta0=[-5.851851	-5.2944001	-2.8233827	-3.196007	-5.7383223	-5.6147304...
 %   -6.8702779	-4.1413453	-4.1940166	-4.7819223	-3.1666896	-3.6302424 ...
  %  1.86681	3.90219	-1.23522  -2.5163	4.01914	4.53675...
  %  4.38982	 -0.996	3.33921	3.28119	2.41559	1.10417...
  %  5.172  2.7550101	-1.49214];   %rescale to reduce rounding error

%theta0=[-5.91740710 -5.37648950 -2.91890680 -3.27645740 -5.79194900 -5.70896380...
  %  -6.77502670 -4.22487880 -4.26353970 -4.90114720 -3.25408260 -3.70270700...
   % 0.00019496 0.00039429 -0.00011557 -0.00023990 0.00038988 0.00045619...
    %0.00037716 -0.00008556 0.00034018 0.00036918 0.00024553 0.00012336...
   % 0.03299912 -1.39623530];
%}
%theta0=[-4.9:-0.1:-6,0.0008:0.0001:0.0019,0.1];
%options = optimset('MaxIter',1e+10000000000000000000000,'MaxFunEvals',1e+10000000000000000000);
options = optimoptions(@fminunc,'MaxIter',1e+10000000000000000000000,'MaxFunEvals',1e+10000000000000000000,'Algorithm','quasi-newton');


%[theta,fval] = fminsearch(@(theta) logL_ReSTAT_mktFEs_pophhi(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,pop65_cub,hhi_cub,df_cub,eerr_cub,pre,act,theta),theta0,options);
[theta,fval] = fminunc(@(theta) logL_ReSTAT(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,pop65_cub,hhi_cub,df_cub,eerr_cub,pre,act,theta),theta0,options);

textFileName_theta=sprintf('theta%d.csv',start);  
dlmwrite(textFileName_theta,[theta fval],'delimiter', ',', 'precision', 8)


%%%%%%%%%%%%%%%%%%
%Compute the standard error
epsilon=1e-6;  %define the disturbance
pr=prob_ReSTAT(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,pop65_cub,hhi_cub,df_cub,eerr_cub,pre,act,theta); %pr26 is the probability function for 26 parameters.
oldpr=pr;
oldlnp=log(pr); %the probability given the original estimates
oldq=theta; %the original estimates
dldq=[];
for i=1:size(theta,2)   %the following for-loop computes the 1st derivative w.r.t. theta
    theta=oldq;
    theta(i)=oldq(i)+epsilon;
    pr1=prob_ReSTAT(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,pop65_cub,hhi_cub,df_cub,eerr_cub,pre,act,theta);
    newlnp=log(pr1);
    dldq=[dldq (newlnp-oldlnp)/epsilon];
end
%A=dldq'*dldq;
%var_q=diag(inv(dldq'*dldq),0);
se_q=sqrt(diag(inv(dldq'*dldq),0));   %the standard error is the sum of the outer product of first derivative

textFileName_se=sprintf('se%d.csv',start);  %the input file  
dlmwrite(textFileName_se,se_q,'delimiter', ',', 'precision', 8)
%}

toc
%diary off

